Preparations

Load the necessary libraries

library(rstanarm)   #for fitting models in STAN
library(brms)       #for fitting models in STAN
library(coda)       #for diagnostics
library(ggmcmc)     #for MCMC diagnostics
library(bayesplot)  #for diagnostics
library(rstan)      #for interfacing with STAN
library(DHARMa)     #for residual diagnostics
library(emmeans)    #for marginal means etc
library(broom)      #for tidying outputs
library(broom.mixed)
library(tidybayes)  #for more tidying outputs
library(ggeffects)  #for partial plots
library(tidyverse)  #for data wrangling etc
library(patchwork)

Scenario

Here is a modified example from Quinn and Keough (2002). Day and Quinn (1989) described an experiment that examined how rock surface type affected the recruitment of barnacles to a rocky shore. The experiment had a single factor, surface type, with 4 treatments or levels: algal species 1 (ALG1), algal species 2 (ALG2), naturally bare surfaces (NB) and artificially scraped bare surfaces (S). There were 5 replicate plots for each surface type and the response (dependent) variable was the number of newly recruited barnacles on each plot after 4 weeks.

Six-plated barnacle

Format of day.csv data files

TREAT	BARNACLE
ALG1	27
..	..
ALG2	24
..	..
NB	9
..	..
S	12
..	..

TREAT	Categorical listing of surface types. ALG1 = algal species 1, ALG2 = algal species 2, NB = naturally bare surface, S = scraped bare surface.
BARNACLE	The number of newly recruited barnacles on each plot after 4 weeks.

Read in the data

day = read_csv('../public/data/day.csv', trim_ws=TRUE)
glimpse(day)

## Rows: 20
## Columns: 2
## $ TREAT    <chr> "ALG1", "ALG1", "ALG1", "ALG1", "ALG1", "ALG2", "ALG2", "ALG…
## $ BARNACLE <dbl> 27, 19, 18, 23, 25, 24, 33, 27, 26, 32, 9, 13, 17, 14, 22, 1…

Start by declaring the categorical variables as factor.

day = day %>% mutate(TREAT=factor(TREAT))

Model formula: \[ \begin{align} y_i &\sim{} \mathcal{Pois}(\lambda_i)\\ ln(\mu_i) &= \boldsymbol{\beta} \bf{X_i}\\ \beta_0 &\sim{} \mathcal{N}(0,10)\\ \beta_{1,2,3} &\sim{} \mathcal{N}(0,1)\\ \end{align} \]

where \(\boldsymbol{\beta}\) is a vector of effects parameters and \(\bf{X}\) is a model matrix representing the intercept and treatment contrasts for the effects of Treatment on barnacle recruitment.

Exploratory data analysis

The exploratory data analyses that we performed in the frequestist installment of this example are equally valid here. That is, boxplots and/or violin plots for each population (substrate type).

ggplot(day, aes(y=BARNACLE, x=TREAT)) +
    geom_boxplot()+
    geom_point(color='red')

ggplot(day, aes(y=BARNACLE, x=TREAT)) +
    geom_violin()+
    geom_point(color='red')

Conclusions:

although exploratory data analysis suggests that we might well be fine modelling these data against a gaussian distribution, a Poisson distribution would clearly be a more natural choice and it would also prevent any non positive predtions.

Fit the model

rstanarm

Using default priors

In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

day.rstanarm = stan_glm(BARNACLE ~ TREAT, data=day,
                        family=poisson(link='log'), 
                         iter = 5000, warmup = 1000,
                         chains = 3, thin = 5, refresh = 0)

prior_summary(day.rstanarm)

## Priors for model 'day.rstanarm' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0], scale = [2.5,2.5,2.5])
##   Adjusted prior:
##     ~ normal(location = [0,0,0], scale = [5.63,5.63,5.63])
## ------
## See help('prior_summary.stanreg') for more details

This tells us:

for the intercept, when the family is Poisson, it is using a normal prior with a mean of 0 and a standard deviation of 2.5. The 2.5 is used for all intercepts. It is often scaled, but only if it is larger than 2.5 is the scaled version kept.
for the coefficients (in this case, just the slope), the default prior is a normal prior centered around 0 with a standard deviation of 2.5. This is then adjusted for the scale of the data by dividing the 2.5 by the standard deviation of the numerical dummy variables for the predictor (then rounded).

2.5/sd(model.matrix(~TREAT, day)[, 2])

## [1] 5.627314

there is no auxillary prior as we are employing a Poisson distribution.

Assessing priors

One way to assess the priors is to have the MCMC sampler sample purely from the prior predictive distribution without conditioning on the observed data. Doing so provides a glimpse at the range of predictions possible under the priors. On the one hand, wide ranging preditions would ensure that the priors are unlikely to influence the actual preditions once they are conditioned on the data. On the other hand, if they are too wide, the sampler is being permitted to traverse into regions of parameter space that are not logically possible in the context of the actual underlying ecological context. Not only could this mean that illogical parameter estimates are possible, when the sampler is traversing regions of paramter space that are not supported by the actual data, the sampler can become unstable and have difficulty.

We can draw from the prior predictive distribution instead of conditioning on the response, by updating the model and indicating prior_PD=TRUE. After refittin the model in this way, we can plot the predictions to gain insights into the range of predictions supported by the priors alone.

day.rstanarm1 <- update(day.rstanarm,  prior_PD=TRUE)

ggpredict(day.rstanarm1) %>% plot(add.data=TRUE)

## $TREAT

Conclusions:

we see that the range of predictions is faily wide and the predicted means could range from 0 to very large (perhaps too large).

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

\(\beta_0\): normal centered at 0 with a standard deviation of 10
\(\beta_1\): normal centered at 0 with a standard deviation of 1

I will also overlay the raw data for comparison.

day.rstanarm2= stan_glm(BARNACLE ~ TREAT, data=day,
                        family=poisson(link='log'), 
                         prior_intercept = normal(0, 10, autoscale=FALSE),
                         prior = normal(0, 1, autoscale=FALSE),
                         prior_PD=TRUE, 
                         iter = 5000, warmup = 1000,
                         chains = 3, thin = 5, refresh = 0
                         )

ggpredict(day.rstanarm2) %>%
  plot(add.data=TRUE)

## $TREAT

Now lets refit, conditioning on the data.

day.rstanarm3= update(day.rstanarm2,  prior_PD=FALSE)

Plotting prior and posterior

posterior_vs_prior(day.rstanarm3, color_by='vs', group_by=TRUE,
                   facet_args=list(scales='free_y'))

Conclusions:

in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

ggpredict(day.rstanarm3) %>% plot(add.data=TRUE)

## $TREAT

brms

Using default priors

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

day.brm = brm(bf(BARNACLE ~ TREAT), data=day,
              family=poisson(link='log'), 
              iter = 5000, warmup = 1000,
              chains = 3, thin = 5, refresh = 0)

prior_summary(day.brm)

##                 prior     class      coef group resp dpar nlpar bound       source
##                (flat)         b                                            default
##                (flat)         b TREATALG2                             (vectorized)
##                (flat)         b   TREATNB                             (vectorized)
##                (flat)         b    TREATS                             (vectorized)
##  student_t(3, 3, 2.5) Intercept                                            default

This tells us:

for the intercept, it is using a student t (flatter normal) prior with a mean of 0 and a standard deviation of 2.5. These mean and standard deviation values are the defaults.
for the beta coefficients (in this case, each effect), the default prior is a inproper flat prior. A flat prior essentially means that any value between negative infinity and positive infinity are equally likely. Whilst this might seem reckless, in practice, it seems to work reasonably well for non-intercept beta parameters.
since we have nominated a Poisson distribution, there is no auxillary prior.

Assessing priors

In brms, we can inform the sampler to draw from the prior predictive distribution instead of conditioning on the response, by running the model with the sample_prior='only' argument. Unfortunately, this cannot be applied when there are flat priors (since the posteriors will necessarily extend to negative and positive infinity). Therefore, in order to use this useful routine, we need to make sure that we have defined a proper prior for all parameters.

day.brm1 = brm(bf(BARNACLE ~ TREAT), data=day,
               family=poisson(link='log'), 
               prior=prior(normal(0, 1), class='b'), 
               sample_prior = 'only', 
               iter = 5000, warmup = 1000,
               chains = 3, thin = 5, refresh = 0)

ggpredict(day.brm1) %>% plot(add.data=TRUE)

## $TREAT

conditional_effects(day.brm1) %>%  plot(points=TRUE)

Conclusions:

we see that the range of predictions is faily wide (ranging from 0 to very high)

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

\(\beta_0\): normal centered at 0 with a standard deviation of 10
\(\beta_1\): normal centered at 0 with a standard deviation of 1

I will also overlay the raw data for comparison.

day.brm2 = brm(bf(BARNACLE ~ TREAT), data=day,
               family=poisson(link='log'), 
               prior=c(prior(normal(0, 10),  class='Intercept'),
                       prior(normal(0, 1), class='b')
                            ), 
                sample_prior = 'only', 
                iter = 5000, warmup = 1000,
                chains = 3, thin = 5, refresh = 0)

ggpredict(day.brm2) %>%
  plot(add.data=TRUE)

## $TREAT

day.brm3 <- update(day.brm2,  sample_prior=TRUE, refresh=0)

Plotting prior and posterior

day.brm3 %>% get_variables()

##  [1] "b_Intercept"     "b_TREATALG2"     "b_TREATNB"       "b_TREATS"       
##  [5] "prior_Intercept" "prior_b"         "lp__"            "accept_stat__"  
##  [9] "stepsize__"      "treedepth__"     "n_leapfrog__"    "divergent__"    
## [13] "energy__"

day.brm3 %>%
  posterior_samples %>%
  select(-`lp__`) %>%
  gather %>%
  mutate(Type=ifelse(str_detect(key, 'prior'), 'Prior', 'b'),
         Class=ifelse(str_detect(key, 'Intercept'),  'Intercept',
               ifelse(str_detect(key, 'sigma'),  'Sigma',  'b'))) %>%
  ggplot(aes(x=Type,  y=value)) +
  stat_pointinterval()+
  facet_wrap(~Class,  scales='free')

Exploring the stan code

standata(day.brm3)

## $N
## [1] 20
## 
## $Y
##  [1] 27 19 18 23 25 24 33 27 26 32  9 13 17 14 22 12  8 15 20 11
## 
## $K
## [1] 4
## 
## $X
##    Intercept TREATALG2 TREATNB TREATS
## 1          1         0       0      0
## 2          1         0       0      0
## 3          1         0       0      0
## 4          1         0       0      0
## 5          1         0       0      0
## 6          1         1       0      0
## 7          1         1       0      0
## 8          1         1       0      0
## 9          1         1       0      0
## 10         1         1       0      0
## 11         1         0       1      0
## 12         1         0       1      0
## 13         1         0       1      0
## 14         1         0       1      0
## 15         1         0       1      0
## 16         1         0       0      1
## 17         1         0       0      1
## 18         1         0       0      1
## 19         1         0       0      1
## 20         1         0       0      1
## attr(,"assign")
## [1] 0 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$TREAT
##      ALG2 NB S
## ALG1    0  0 0
## ALG2    1  0 0
## NB      0  1 0
## S       0  0 1
## 
## 
## $prior_only
## [1] 0
## 
## attr(,"class")
## [1] "standata"

stancode(day.brm3)

## // generated with brms 2.14.0
## functions {
## }
## data {
##   int<lower=1> N;  // total number of observations
##   int Y[N];  // response variable
##   int<lower=1> K;  // number of population-level effects
##   matrix[N, K] X;  // population-level design matrix
##   int prior_only;  // should the likelihood be ignored?
## }
## transformed data {
##   int Kc = K - 1;
##   matrix[N, Kc] Xc;  // centered version of X without an intercept
##   vector[Kc] means_X;  // column means of X before centering
##   for (i in 2:K) {
##     means_X[i - 1] = mean(X[, i]);
##     Xc[, i - 1] = X[, i] - means_X[i - 1];
##   }
## }
## parameters {
##   vector[Kc] b;  // population-level effects
##   real Intercept;  // temporary intercept for centered predictors
## }
## transformed parameters {
## }
## model {
##   // likelihood including all constants
##   if (!prior_only) {
##     target += poisson_log_glm_lpmf(Y | Xc, Intercept, b);
##   }
##   // priors including all constants
##   target += normal_lpdf(b | 0, 1);
##   target += normal_lpdf(Intercept | 0, 10);
## }
## generated quantities {
##   // actual population-level intercept
##   real b_Intercept = Intercept - dot_product(means_X, b);
##   // additionally draw samples from priors
##   real prior_b = normal_rng(0,1);
##   real prior_Intercept = normal_rng(0,10);
## }

MCMC sampling diagnostics

In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.

There are a wide variety of tests that range from the big picture, overal chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.

rstanarm

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

available_mcmc()

## bayesplot MCMC module:
##   mcmc_acf
##   mcmc_acf_bar
##   mcmc_areas
##   mcmc_areas_data
##   mcmc_areas_ridges
##   mcmc_areas_ridges_data
##   mcmc_combo
##   mcmc_dens
##   mcmc_dens_chains
##   mcmc_dens_chains_data
##   mcmc_dens_overlay
##   mcmc_hex
##   mcmc_hist
##   mcmc_hist_by_chain
##   mcmc_intervals
##   mcmc_intervals_data
##   mcmc_neff
##   mcmc_neff_data
##   mcmc_neff_hist
##   mcmc_nuts_acceptance
##   mcmc_nuts_divergence
##   mcmc_nuts_energy
##   mcmc_nuts_stepsize
##   mcmc_nuts_treedepth
##   mcmc_pairs
##   mcmc_parcoord
##   mcmc_parcoord_data
##   mcmc_rank_hist
##   mcmc_rank_overlay
##   mcmc_recover_hist
##   mcmc_recover_intervals
##   mcmc_recover_scatter
##   mcmc_rhat
##   mcmc_rhat_data
##   mcmc_rhat_hist
##   mcmc_scatter
##   mcmc_trace
##   mcmc_trace_data
##   mcmc_trace_highlight
##   mcmc_violin

Of these, we will focus on:

mcmc_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

plot(day.rstanarm3, plotfun='mcmc_trace')

The chains appear well mixed and very similar

acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

plot(day.rstanarm3, 'acf_bar')

There is no evidence of autocorrelation in the MCMC samples

Rhat: Rhat is a measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

plot(day.rstanarm3, 'rhat_hist')

All Rhat values are below 1.05, suggesting the chains have converged.

neff (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

plot(day.rstanarm3, 'neff_hist')

Ratios all very high.

More diagnostics

plot(day.rstanarm3, 'combo')

plot(day.rstanarm3, 'violin')

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

stan_trace(day.rstanarm3)

The chains appear well mixed and very similar

stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

stan_ac(day.rstanarm3)

There is no evidence of autocorrelation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

stan_rhat(day.rstanarm3)

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

stan_ess(day.rstanarm3)

Ratios all very high.

stan_dens(day.rstanarm3, separate_chains = TRUE)

ggmcmc

The ggmean package also has a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

day.ggs <- ggs(day.rstanarm3)
ggs_traceplot(day.ggs)

The chains appear well mixed and very similar

gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

ggs_autocorrelation(day.ggs)

There is no evidence of autocorrelation in the MCMC samples

stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

ggs_Rhat(day.ggs)

All Rhat values are below 1.05, suggesting the chains have converged.

stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

If the ratios are low, tightening the priors may help.

ggs_effective(day.ggs)

Ratios all very high.

More diagnostics

ggs_crosscorrelation(day.ggs)

ggs_grb(day.ggs)

Model validation

rstanarm

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

available_ppc()

## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_data
##   ppc_dens
##   ppc_dens_overlay
##   ppc_ecdf_overlay
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_data
##   ppc_intervals_grouped
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_data
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped

dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

pp_check(day.rstanarm3,  plotfun='dens_overlay')

The model draws appear deviate from the observed data.

error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.

pp_check(day.rstanarm3, plotfun='error_scatter_avg')

The predictive error seems to be related to the predictor - the model performs poorest at higher mussel clump areas.

error_scatter_avg_vs_x: this is similar to a regular residual plot and as such should be interpreted as such. Again, this is not interpretable for binary data.

pp_check(day.rstanarm3, x=as.numeric(day$TREAT), plotfun='error_scatter_avg_vs_x')

intervals: plots the observed data overlayed ontop of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

pp_check(day.rstanarm3, x=as.numeric(day$TREAT), plotfun='intervals')

The modelled preditions seem to underestimate the uncertainty with increasing mussel clump area.

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

#library(shinystan)
#launch_shinystan(day.rstanarm3)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components ourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

simulated (predicted) responses associated with each observation.
observed values
fitted (predicted) responses (averaged) associated with each observation

preds <- posterior_predict(day.rstanarm3,  nsamples=250,  summary=FALSE)
day.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = day$BARNACLE,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse = TRUE)
plot(day.resids)

Conclusions:

the simulated residuals DOES NOT suggest any issues with the fitted model

brms

pp check

available_ppc()

## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_data
##   ppc_dens
##   ppc_dens_overlay
##   ppc_ecdf_overlay
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_data
##   ppc_intervals_grouped
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_data
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped

dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

pp_check(day.brm3,  type='dens_overlay')

The model draws appear deviate from the observed data.

error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. There is some pattern remaining in these residuals.

pp_check(day.brm3, type='error_scatter_avg')

The predictive error seems to be related to the predictor - the model performs poorest at higher mussel clump areas.

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

#library(shinystan)
#launch_shinystan(day.brm3)

DHARMa residuals

We need to supply:

simulated (predicted) responses associated with each observation.
observed values
fitted (predicted) responses (averaged) associated with each observation

preds <- posterior_predict(day.brm3,  nsamples=250,  summary=FALSE)
day.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = day$BARNACLE,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse = TRUE)
plot(day.resids)

Conclusions:

the simulated residuals DO NOT suggest any issues with the model fit

Partial effects plots

rstanarm

ggpredict

day.rstanarm3 %>% ggpredict() %>% plot(add.data=TRUE)

## $TREAT

ggemmeans

day.rstanarm3 %>% ggemmeans(~TREAT,  type='fixed') %>% plot(add.data=TRUE)

fitted_draws

day.rstanarm3 %>% fitted_draws(newdata=day) %>%
  median_hdci() %>%
  ggplot(aes(x=TREAT, y=.value)) +
  geom_pointrange(aes(ymin=.lower, ymax=.upper)) + 
  geom_line() +
  geom_point(data=day,  aes(y=BARNACLE,  x=TREAT))

brms

conditional_effects

day.brm3 %>% conditional_effects()

ggpredict

day.brm3 %>% ggpredict() %>% plot(add.data=TRUE)

## $TREAT

ggemmeans

day.brm3 %>% ggemmeans(~TREAT) %>% plot(add.data=TRUE)

fitted_draws

day.brm3 %>% fitted_draws(newdata=day) %>%
  median_hdci() %>%
  ggplot(aes(x=TREAT, y=.value)) +
  geom_pointrange(aes(ymin=.lower, ymax=.upper)) + 
  geom_line() +
  geom_point(data=day,  aes(y=BARNACLE,  x=TREAT))

Model investigation

rstanarm

rstanarm captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

summary(day.rstanarm3)

## 
## Model Info:
##  function:     stan_glm
##  family:       poisson [log]
##  formula:      BARNACLE ~ TREAT
##  algorithm:    sampling
##  sample:       2400 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 20
##  predictors:   4
## 
## Estimates:
##               mean   sd   10%   50%   90%
## (Intercept)  3.1    0.1  3.0   3.1   3.2 
## TREATALG2    0.2    0.1  0.1   0.2   0.4 
## TREATNB     -0.4    0.1 -0.6  -0.4  -0.2 
## TREATS      -0.5    0.2 -0.7  -0.5  -0.3 
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 19.8    1.4 18.0  19.8  21.6 
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.0  1.0  2438 
## TREATALG2     0.0  1.0  2312 
## TREATNB       0.0  1.0  2425 
## TREATS        0.0  1.0  2518 
## mean_PPD      0.0  1.0  2202 
## log-posterior 0.0  1.0  2301 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Conclusions:

in the Model Info, we are informed that the total MCMC posterior sample size is 2400 and that there were 20 raw observations.
the estimated mean (expected number of newly recruited barnacles) on the ALG1 surface is 3.1. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 22.15.
the estimated effect of ALG2 vs ALG1 is 0.24 (mean) or 0.09 (median) with a standard deviation of 0. The 90% credibility intervals indicate that we are 90% confident that the slope is between 0.24 and 0.24 - e.g. there is a significant positive effect. When back-transformed onto the reponse scale, we see that barnacle recruitment on ALG2 is 1.27 times higher than that on ALG1. This represents a 27% increase in barnacle recruitment.
the estimated effect of NB and S are -0.39 and -0.52 respectively, which equate to 1.48 and 1.68 fold declines respectively.
Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

tidyMCMC

tidyMCMC(day.rstanarm3$stanfit, estimate.method='median',  conf.int=TRUE,  conf.method='HPDinterval',  rhat=TRUE, ess=TRUE)

Conclusions:

the estimated mean (expected number of newly recruited barnacles) on the ALG1 surface is 3.1. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 22.19.
the estimated effect of ALG2 vs ALG1 is 0.24 (median) with a standard error of 0.12. The 95% credibility intervals indicate that we are 95% confident that the effect is between 0.01 and 0.49 - e.g. there is a significant positive effect. When back-transformed onto the reponse scale, we see that barnacle recruitment on ALG2 is 1.28 times higher than that on ALG1. This represents a 28% increase in barnacle recruitment.
the estimated effect of NB and S are -0.39 and -0.52 respectively, which equate to 1.48 and 1.68 fold declines respectively.
Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

gather_draws

Due to the presence of a log transform in the predictor, it is better to use the regex version.

day.rstanarm3 %>% get_variables()

##  [1] "(Intercept)"   "TREATALG2"     "TREATNB"       "TREATS"       
##  [5] "accept_stat__" "stepsize__"    "treedepth__"   "n_leapfrog__" 
##  [9] "divergent__"   "energy__"

day.draw <- day.rstanarm3 %>% gather_draws(`.Intercept.*|.*TREAT.*`,  regex=TRUE)
day.draw

We can then summarise this

day.draw %>% median_hdci

day.rstanarm3 %>% 
  gather_draws(`.Intercept.*|.*TREAT.*`, regex=TRUE) %>% 
  ggplot() + 
  stat_halfeye(aes(x=.value,  y=.variable)) +
  facet_wrap(~.variable, scales='free')

We could alternatively express the parameters on the response scale.

day.rstanarm3 %>% 
  gather_draws(`.Intercept.*|.*TREAT.*`, regex=TRUE) %>%
  group_by(.variable) %>%
  mutate(.value=exp(.value)) %>%
  median_hdci

Conclusions:

the estimated mean (expected number of newly recruited barnacles) on the ALG1 surface is 3.1. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 22.19.
the estimated effect of ALG2 vs ALG1 is 0.24 (median) with a standard error of 0.01. The 95% credibility intervals indicate that we are 95% confident that the effect is between 0.49 and 0.95 - e.g. there is a significant positive effect. When back-transformed onto the reponse scale, we see that barnacle recruitment on ALG2 is 1.28 times higher than that on ALG1. This represents a 28% increase in barnacle recruitment.
the estimated effect of NB and S are -0.39 and -0.52 respectively, which equate to 1.48 and 1.68 fold declines respectively.
Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

bayesplot

day.rstanarm3 %>% plot(plotfun='mcmc_intervals')

tidy_draws

This is purely a graphical depiction on the posteriors.

day.rstanarm3 %>% tidy_draws()

spread_draws

day.rstanarm3 %>% spread_draws(`.Intercept.*|.*TREAT.*`,  regex=TRUE)

posterior_samples

day.rstanarm3 %>% posterior_samples() %>% as_tibble()

\(R^2\)

Unfortunately, \(R^2\) calculations for models other than Gaussian and Binomial have not yet been implemented for rstanarm models yet.

#day.rstanarm3 %>% bayes_R2() %>% median_hdci

brms

brms captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

summary(day.brm3)

##  Family: poisson 
##   Links: mu = log 
## Formula: BARNACLE ~ TREAT 
##    Data: day (Number of observations: 20) 
## Samples: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup samples = 2400
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     3.10      0.09     2.91     3.28 1.00     2297     2053
## TREATALG2     0.25      0.13    -0.01     0.50 1.00     2075     1898
## TREATNB      -0.39      0.15    -0.69    -0.11 1.00     2346     2329
## TREATS       -0.52      0.16    -0.84    -0.22 1.00     2306     2082
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Conclusions:

in the Model Info, we are informed that the total MCMC posterior sample size is 2400 and that there were 20 raw observations.
the estimated mean (expected number of newly recruited barnacles) on the ALG1 surface is 3.1. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 22.14.
the estimated effect of ALG2 vs ALG1 is 0.25 (mean) or 0.5 (median) with a standard deviation of 0.13. The 90% credibility intervals indicate that we are 90% confident that the slope is between 0.25 and 1 - e.g. there is a significant positive effect. When back-transformed onto the reponse scale, we see that barnacle recruitment on ALG2 is 1.28 times higher than that on ALG1. This represents a 28% increase in barnacle recruitment.
the estimated effect of NB and S are -0.39 and -0.52 respectively, which equate to 1.48 and 1.68 fold declines respectively.
Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

tidyMCMC

tidyMCMC(day.brm3$fit, estimate.method='median',  conf.int=TRUE,  conf.method='HPDinterval',  rhat=TRUE, ess=TRUE)

Conclusions:

the estimated mean (expected number of newly recruited barnacles) on the ALG1 surface is 3.1. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 22.14.
the estimated effect of ALG2 vs ALG1 is 0.25 (median) with a standard error of 0.13. The 95% credibility intervals indicate that we are 95% confident that the effect is between -0.01 and 0.5 - e.g. there is a significant positive effect. When back-transformed onto the reponse scale, we see that barnacle recruitment on ALG2 is 1.28 times higher than that on ALG1. This represents a 28% increase in barnacle recruitment.
the estimated effect of NB and S are -0.39 and -0.51 respectively, which equate to 1.47 and 1.67 fold declines respectively.
Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

gather_draws

Due to the presence of a log transform in the predictor, it is better to use the regex version.

day.brm3 %>% get_variables()

##  [1] "b_Intercept"     "b_TREATALG2"     "b_TREATNB"       "b_TREATS"       
##  [5] "prior_Intercept" "prior_b"         "lp__"            "accept_stat__"  
##  [9] "stepsize__"      "treedepth__"     "n_leapfrog__"    "divergent__"    
## [13] "energy__"

day.draw <- day.brm3 %>% gather_draws(`.Intercept.*|.*TREAT.*`,  regex=TRUE)
day.draw

We can then summarise this

day.draw %>% median_hdci

day.brm3 %>% 
  gather_draws(`.Intercept.*|.*TREAT.*`, regex=TRUE) %>% 
  ggplot() + 
  stat_halfeye(aes(x=.value,  y=.variable)) +
  facet_wrap(~.variable, scales='free')

We could alternatively express the parameters on the response scale.

day.brm3 %>% 
  gather_draws(`.Intercept.*|.*TREAT.*`, regex=TRUE) %>%
  group_by(.variable) %>%
  mutate(.value=exp(.value)) %>%
  median_hdci

Conclusions:

the estimated mean (expected number of newly recruited barnacles) on the ALG1 surface is 0.25. This is the mean of the posterior distribution for this parameter. If we back-transform this to the response scale, this becomes 1.28.
the estimated effect of ALG2 vs ALG1 is -0.39 (median) with a standard error of -0.69. The 95% credibility intervals indicate that we are 95% confident that the effect is between -0.11 and 0.95 - e.g. there is a significant positive effect. When back-transformed onto the reponse scale, we see that barnacle recruitment on ALG2 is 0.68 times higher than that on ALG1. This represents a -32% increase in barnacle recruitment.
the estimated effect of NB and S are -0.51 and NA respectively, which equate to 1.67 and NA fold declines respectively.
Rhat and number of effective samples for each parameter are also provided as MCMC diagnostics and all look good.

bayesplot

day.brm3 %>% plot(plotfun='mcmc_intervals')

tidy_draws

This is purely a graphical depiction on the posteriors.

day.brm3 %>% tidy_draws()

spread_draws

day.brm3 %>% spread_draws(`.Intercept.*|.*TREAT.*`,  regex=TRUE)

posterior_samples

day.brm3 %>% posterior_samples() %>% as_tibble()

\(R^2\)

day.brm3 %>% bayes_R2(summary=FALSE) %>% median_hdci

Further investigations

The estimated coefficients as presented in the summary tables above highlight very specific comparisons. However, there are other possible comparisons that we might be interested in. For example, whist the treatment effects compare each of the substrate types against the first level (ALG1), we might also be interested in the differences between (for example) the two bare substrates (NB and S).

To get at more of the comparisons we have two broad approaches:

compare all substrates to each other in a pairwise manner.
- importantly, given that the probability of a Type I error (falsely rejecting a null hypothesis) is set at 0.05 per comparison, when there are multiple comparisons, the family-wise Type I error (probability of at least one false rejection amongst all comparisons) rate compounds. It is important to attetmpt to constrain this family-wise Type I error rate. One approach to do this is a Tukey’s test.
- keep in mind that in order to constrain the family-wise Type I error rate at 0.05, the power of each individual comparison is reduced (individual Type II error rates increased).
define a small set of specific comparisons. There should be a maximum of \(p-1\) comparisons defined (where \(p\) is the number of levels of the categorical predictors) and each comparison should be independent of one another.

We will now only persue the Poisson model.

Post-hoc test (Tukey’s)

rstanarm

effects

day.rstanarm3 %>% emmeans(pairwise ~TREAT, type='response')

## $emmeans
##  TREAT rate lower.HPD upper.HPD
##  ALG1  22.2      18.4      26.3
##  ALG2  28.2      23.6      32.8
##  NB    15.0      11.9      18.3
##  S     13.3      10.2      16.8
## 
## Point estimate displayed: median 
## Results are back-transformed from the log scale 
## HPD interval probability: 0.95 
## 
## $contrasts
##  contrast    ratio lower.HPD upper.HPD
##  ALG1 / ALG2 0.783     0.601     0.973
##  ALG1 / NB   1.480     1.091     1.921
##  ALG1 / S    1.675     1.187     2.220
##  ALG2 / NB   1.886     1.416     2.421
##  ALG2 / S    2.129     1.559     2.822
##  NB / S      1.135     0.782     1.543
## 
## Point estimate displayed: median 
## Results are back-transformed from the log scale 
## HPD interval probability: 0.95

Conclusions:

the contrasts section of the output indicates that there is evidence of:
ALG1 has 1.48 fold (48%) more newly recruited barnacles than the NB substrate
ALG1 has 1.68 fold (68%) more newly recruited barnacles than the S substrate
ALG2 has 1.89 fold (89%) more newly recruited barnacles than the NB substrate
ALG2 has 2.13 fold (113%) more newly recruited barnacles than the S substrate
ALG1 was not found to be different to ALG2 and NB was not found to be different to S

Probabilities

day.em = emmeans(day.rstanarm3, pairwise~TREAT, type='link')$contrasts %>%
    gather_emmeans_draws() %>%
    mutate(Fit=exp(.value))
day.em %>% head

day.em %>% group_by(contrast) %>%
    ggplot(aes(x=Fit)) +
    geom_histogram() +
    geom_vline(xintercept=1, color='red') + 
    facet_wrap(~contrast, scales='free')

day.em %>% group_by(contrast) %>% median_hdi()

# Probability of effect
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1)/n())

##Probability of effect greater than 10%
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1.1)/n())

brms

effects

day.brm3 %>% emmeans(pairwise~TREAT, type='response')

## $emmeans
##  TREAT rate lower.HPD upper.HPD
##  ALG1  22.1      18.1      26.2
##  ALG2  28.3      23.6      33.2
##  NB    15.0      11.8      18.4
##  S     13.3      10.2      16.7
## 
## Point estimate displayed: median 
## Results are back-transformed from the log scale 
## HPD interval probability: 0.95 
## 
## $contrasts
##  contrast    ratio lower.HPD upper.HPD
##  ALG1 / ALG2 0.783     0.603      1.00
##  ALG1 / NB   1.474     1.078      1.96
##  ALG1 / S    1.672     1.221      2.25
##  ALG2 / NB   1.885     1.370      2.44
##  ALG2 / S    2.133     1.539      2.81
##  NB / S      1.136     0.796      1.54
## 
## Point estimate displayed: median 
## Results are back-transformed from the log scale 
## HPD interval probability: 0.95

Conclusions:

the contrasts section of the output indicates that there is evidence of:
ALG1 has 1.47 fold (47%) more newly recruited barnacles than the NB substrate
ALG1 has 1.67 fold (67%) more newly recruited barnacles than the S substrate
ALG2 has 1.89 fold (89%) more newly recruited barnacles than the NB substrate
ALG2 has 2.13 fold (113%) more newly recruited barnacles than the S substrate
ALG1 was not found to be different to ALG2 and NB was not found to be different to S

Probabilities

day.em = emmeans(day.brm3, pairwise~TREAT, type='link')$contrasts %>%
    gather_emmeans_draws() %>%
    mutate(Fit=exp(.value))
day.em %>% head

day.em %>% group_by(contrast) %>%
    ggplot(aes(x=Fit)) +
    geom_histogram() +
    geom_vline(xintercept=1, color='red') + 
    facet_wrap(~contrast, scales='free')

day.em %>% group_by(contrast) %>% median_hdi()

# Probability of effect
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1)/n())

##Probability of effect greater than 10%
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1.1)/n())

Planned contrasts

Define your own

Compare:

ALG1 vs ALG2
NB vs S
average of ALG1+ALG2 vs NB+S

Levels	Alg1 vs Alg2	NB vs S	Alg vs Bare
Alg1	1	0	0.5
Alg2	-1	0	0.5
NB	0	1	-0.5
S	0	-1	-0.5

rstanarm

##Planned contrasts

## $emmeans
##  TREAT emmean lower.HPD upper.HPD
##  ALG1    3.10      2.92      3.28
##  ALG2    3.34      3.18      3.51
##  NB      2.71      2.49      2.93
##  S       2.59      2.33      2.83
## 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.95 
## 
## $contrasts
##  contrast        estimate lower.HPD upper.HPD
##  TREAT.Alg2_Alg1    0.244   0.00505     0.486
##  TREAT.NB_S         0.127  -0.18410     0.489
##  TREAT.Alg_Bare     0.576   0.38315     0.784
##  TREAT.Alg_NB       0.513   0.24540     0.745
## 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.95

## $emmeans
##  TREAT rate lower.HPD upper.HPD
##  ALG1  22.2      18.4      26.3
##  ALG2  28.2      23.6      32.8
##  NB    15.0      11.9      18.3
##  S     13.3      10.2      16.8
## 
## Point estimate displayed: median 
## Results are back-transformed from the log scale 
## HPD interval probability: 0.95 
## 
## $contrasts
##  contrast        ratio lower.HPD upper.HPD
##  TREAT.Alg2_Alg1  1.28     0.994      1.61
##  TREAT.NB_S       1.13     0.782      1.54
##  TREAT.Alg_Bare   1.78     1.439      2.16
##  TREAT.Alg_NB     1.67     1.278      2.11
## 
## Point estimate displayed: median 
## Results are back-transformed from the log scale 
## HPD interval probability: 0.95

brms

##Planned contrasts

## $emmeans
##  TREAT emmean lower.HPD upper.HPD
##  ALG1    3.10      2.90      3.27
##  ALG2    3.34      3.18      3.51
##  NB      2.71      2.47      2.92
##  S       2.58      2.34      2.83
## 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.95 
## 
## $contrasts
##  contrast        estimate lower.HPD upper.HPD
##  TREAT.Alg2_Alg1    0.245  -0.00744     0.497
##  TREAT.NB_S         0.128  -0.20155     0.457
##  TREAT.Alg_Bare     0.575   0.37553     0.782
##  TREAT.Alg_NB       0.510   0.27110     0.778
## 
## Point estimate displayed: median 
## Results are given on the log (not the response) scale. 
## HPD interval probability: 0.95

## $emmeans
##  TREAT rate lower.HPD upper.HPD
##  ALG1  22.1      18.1      26.2
##  ALG2  28.3      23.6      33.2
##  NB    15.0      11.8      18.4
##  S     13.3      10.2      16.7
## 
## Point estimate displayed: median 
## Results are back-transformed from the log scale 
## HPD interval probability: 0.95 
## 
## $contrasts
##  contrast        ratio lower.HPD upper.HPD
##  TREAT.Alg2_Alg1  1.28     0.965      1.61
##  TREAT.NB_S       1.14     0.796      1.54
##  TREAT.Alg_Bare   1.78     1.447      2.17
##  TREAT.Alg_NB     1.67     1.311      2.18
## 
## Point estimate displayed: median 
## Results are back-transformed from the log scale 
## HPD interval probability: 0.95

Summary Figure

rstanarm

newdata = emmeans(day.rstanarm3, ~TREAT, type='response') %>%
    as.data.frame
newdata

## A quick version
g2 <- ggplot(newdata, aes(y=rate, x=TREAT)) +
    geom_pointrange(aes(ymin=lower.HPD, ymax=upper.HPD)) +
    theme_classic()
    
g2 + g1

brms

newdata = emmeans(day.brm3, ~TREAT, type='response') %>%
    as.data.frame
newdata

## A quick version
g2 <- ggplot(newdata, aes(y=rate, x=TREAT)) +
    geom_pointrange(aes(ymin=lower.HPD, ymax=upper.HPD)) +
    theme_classic()
    
g2 + g1

day.rstanarm <- stan_glm(BARNACLE ~ TREAT, data=day,
                      family=poisson(link='log'),
                      chains = 3,iter = 5000, warmup=2000, thin=5,
                      refresh=0)
prior_summary(day.rstanarm)

day.rstanarm <- stan_glm(BARNACLE ~ TREAT, data=day, family='poisson',
                      prior = normal(c(0,0,0), c(2.5,2.5,2.5)),
                      prior_intercept = normal(0,10),
                      chains = 3,iter = 2000, thin=2, refresh=0)
prior_summary(day.rstanarm)

plot(day.rstanarm,  'mcmc_trace')
plot(day.rstanarm,  'mcmc_acf_bar')
plot(day.rstanarm,  'mcmc_rhat_hist')
plot(day.rstanarm,  'mcmc_neff_hist')


preds <- posterior_predict(day.rstanarm,  nsamples=250,  summary=FALSE)
day.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = day$BARNACLE,
                           fittedPredictedResponse = apply(preds, 2, median),
                           integerResponse = TRUE)
plot(day.resids)


#pp_check(day.rstanarm, x=as.numeric(day$TREAT),'intervals')




## Compare the proportion of zeros in the observed and expected data
#yrep = posterior_predict(day.rstan)
prop_zero <- function(y) mean(y == 0)
(prop_zero_test1 <- pp_check(day.rstanarm, plotfun = "stat", stat = "prop_zero"))
                                        # no zeros - so not zero inflated


day.rstanarmNB <- stan_glm(BARNACLE ~ TREAT, data=day,
                      family='neg_binomial_2',
                      chains = 3,iter = 5000, thin=5, warmup=2000, refresh=0)
preds <- posterior_predict(day.rstanarmNB,  nsamples=250,  summary=FALSE)
day.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = day$BARNACLE,
                            fittedPredictedResponse = apply(preds, 2, median))
plot(day.resids)


day.rstanNB <- update(day.rstan, family = neg_binomial_2)
(loo.P=loo(day.rstan))
(loo.NB=loo(day.rstanNB))
compare_models(loo.P, loo.NB)

ggpredict(day.rstanarm, term='TREAT') %>% plot
ggpredict(day.rstanarm, ~TREAT) %>% plot
ggemmeans(day.rstanarm, ~TREAT) %>% plot


summary(day.rstanarm)
library(tidybayes)
tidyMCMC(day.rstanarm$stanfit, conf.int=TRUE,
         conf.method='HPDinterval', rhat=TRUE,ess=TRUE)


# Pairwise comparisons
library(emmeans)
## factor statements
emmeans(day.rstanarm, pairwise~TREAT, type='response')
## what about probabilities
day.em = emmeans(day.rstanarm, pairwise~TREAT, type='link')$contrasts %>%
    gather_emmeans_draws() %>%
    mutate(Fit=exp(.value))
day.em %>% head
day.em %>% group_by(contrast) %>%
    ggplot(aes(x=Fit)) +
    geom_histogram() +
    geom_vline(xintercept=1, color='red') + 
    facet_wrap(~contrast, scales='free')
day.em %>% group_by(contrast) %>% median_hdi()
# Probability of effect
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1)/n())
##Probability of effect greater than 10%
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1.1)/n())


##Planned contrasts
cmat<-cbind('Alg2_Alg1'=c(-1,1,0,0),
              'NB_S'=c(0,0,1,-1),
             'Alg_Bare'=c(0.5,0.5,-0.5,-0.5),
             'Alg_NB'=c(0.5,0.5,-1,0))
#crossprod(cmat)
emmeans(day.rstanarm, ~TREAT, contr=list(TREAT=cmat), type='link')
emmeans(day.rstanarm, ~TREAT, contr=list(TREAT=cmat), type='response')
day.em = emmeans(day.rstanarm, ~TREAT, contr=list(TREAT=cmat), type='link')$contrasts %>%
      gather_emmeans_draws() %>% mutate(Fit=exp(.value)) 
day.em %>% group_by(contrast) %>% mean_hdi()
# Probability of effect
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1)/n())
##Probability of effect greater than 10%
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1.5)/n())

hist(bayes_R2(day.rstanarmNB))

bayes_R2(day.rstanarm) %>% median_hdi
bayes_R2(day.rstanarmNB) %>% hist


## Summary plot
day.grid = with(day, list(TREAT=levels(TREAT)))
newdata = emmeans(day.rstanarm, ~TREAT, type='response') %>% as.data.frame
head(newdata)
ggplot(newdata, aes(y=rate, x=TREAT)) +
    geom_pointrange(aes(ymin=lower.HPD, ymax=upper.HPD))

day.form <- bf(BARNACLE ~ TREAT,  family=poisson(link='log'))
get_prior(day.form,  data=day)
day.priors <- c(
  prior(normal(0, 10),  class='Intercept'),
  prior(normal(0, 2.5), class='b')
)
day.brms <- brm(day.form, data=day,
                prior=day.priors, 
                 chains=3,  iter=5000,  warmup=2000, thin=5,
                 refresh=0)

plot(day.brms)
mcmc_plot(day.brms,  type='acf_bar')
mcmc_plot(day.brms,  type='rhat_hist')
mcmc_plot(day.brms,  type='neff_hist')

preds <- posterior_predict(day.brms,  nsamples=250,  summary=FALSE)
day.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = day$BARNACLE,
                           fittedPredictedResponse = apply(preds, 2, median),
                           integerResponse = TRUE)
plot(day.resids)

                                        #pp_check(day.brmsP, x=as.numeric(day$TREAT),'intervals')


ggpredict(day.brms, term='TREAT') %>% plot
ggpredict(day.brms, ~TREAT) %>% plot
ggemmeans(day.brms, ~TREAT) %>% plot

summary(day.brms)

tidyMCMC(day.brms$fit, conf.int=TRUE,
         conf.method='HPDinterval', rhat=TRUE,ess=TRUE)

# Pairwise comparisons
library(emmeans)
## factor statements
emmeans(day.brms, pairwise~TREAT, type='response')
## what about probabilities
day.em = emmeans(day.brms, pairwise~TREAT, type='link')$contrasts %>%
    gather_emmeans_draws() %>%
    mutate(Fit=exp(.value))
day.em %>% head
day.em %>% group_by(contrast) %>%
    ggplot(aes(x=Fit)) +
    geom_histogram() +
    geom_vline(xintercept=1, color='red') + 
    facet_wrap(~contrast, scales='free')
day.em %>% group_by(contrast) %>% median_hdi(.width=c(0.8, 0.95))

day.sum <- day.em %>%
  group_by(contrast) %>%
  median_hdci(.width=c(0.8, 0.95))
day.sum
ggplot(day.sum) +
  geom_hline(yintercept=1, linetype='dashed') +
  geom_pointrange(aes(x=contrast, y=Fit, ymin=Fit.lower, ymax=Fit.upper, size=factor(.width)),
                  show.legend = FALSE) +
  scale_size_manual(values=c(1, 0.5)) +
  coord_flip()

g1 <- ggplot(day.sum) +
  geom_hline(yintercept=1) +
  geom_pointrange(aes(x=contrast, y=Fit, ymin=Fit.lower, ymax=Fit.upper, size=factor(.width)), show.legend = FALSE) +
  scale_size_manual(values=c(1, 0.5)) +
  scale_y_continuous(trans=scales::log2_trans(),  breaks=c(0.5, 1, 2, 4)) +
  coord_flip()
g1
                                        # Probability of effect
day.em %>% group_by(contrast) %>% summarize(P=sum(.value>0)/n())
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1)/n())
##Probability of effect greater than 10%
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1.1)/n())

##Planned contrasts
cmat<-cbind('Alg2_Alg1'=c(-1,1,0,0),
              'NB_S'=c(0,0,1,-1),
             'Alg_Bare'=c(0.5,0.5,-0.5,-0.5),
             'Alg_NB'=c(0.5,0.5,-1,0))
#crossprod(cmat)
emmeans(day.brms, ~TREAT, contr=list(TREAT=cmat), type='link')
emmeans(day.brms, ~TREAT, contr=list(TREAT=cmat), type='response')
day.em = emmeans(day.brms, ~TREAT, contr=list(TREAT=cmat), type='link')$contrasts %>%
      gather_emmeans_draws() %>%
      mutate(Fit=exp(.value)) 
day.em %>% group_by(contrast) %>% median_hdci()
# Probability of effect
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1)/n())
##Probability of effect greater than 10%
day.em %>% group_by(contrast) %>% summarize(P=sum(Fit>1.1)/n())

hist(bayes_R2(day.brms, summary=FALSE))

bayes_R2(day.brms, summary=FALSE) %>% median_hdi


## Summary plot
day.grid = with(day, list(TREAT=levels(TREAT)))
newdata = emmeans(day.brms, ~TREAT, type='response') %>% as.data.frame
head(newdata)
g2 <- ggplot(newdata, aes(y=rate, x=TREAT)) +
  geom_pointrange(aes(ymin=lower.HPD, ymax=upper.HPD))

library(patchwork)
g1 + g2

References

Quinn, G. P., and K. J. Keough. 2002. Experimental Design and Data Analysis for Biologists. London: Cambridge University Press.

Bayesian GLM Part5

Murray Logan

12 December, 2020

Preparations

Scenario

Read in the data

Exploratory data analysis

Fit the model

rstanarm

Using default priors

Assessing priors

Defining priors

Plotting prior and posterior

brms

Using default priors

Assessing priors

Defining priors

Plotting prior and posterior

Exploring the stan code

MCMC sampling diagnostics

rstanarm

bayesplot

stan plots

ggmcmc

Model validation

rstanarm

pp check

DHARMa residuals

brms

pp check

DHARMa residuals

Partial effects plots

rstanarm

ggpredict

ggemmeans

fitted_draws

brms

conditional_effects

ggpredict

ggemmeans

fitted_draws

Model investigation

rstanarm

summary

tidyMCMC

gather_draws

bayesplot

tidy_draws

spread_draws

posterior_samples

\(R^2\)

brms

summary

tidyMCMC

gather_draws

bayesplot

tidy_draws

spread_draws

posterior_samples

\(R^2\)

Further investigations

Post-hoc test (Tukey’s)

rstanarm

effects

Probabilities

brms

effects

Probabilities

Planned contrasts

rstanarm

brms

Summary Figure

rstanarm

brms

References